MatLab等值面函数的Python/Numpy等价

2023-03-15 17:30| 来源: 网络整理| 查看: 265

尽管PyVista不在您的目标库中，但基于VTK构建的PyVista可以帮助您轻松完成这项工作(免责声明：我是开发人员之一)。既然你在评论中似乎接受了基于PyVista的解决方案，下面是你会怎么做：1.为您的数据类型定义网格，通常为StructuredGrid，尽管示例中的等距网格甚至可以与UniformGrid一起使用，1.用contour滤镜计算其等值面，1.使用包含等值面的网格的save方法另存为.stl文件。

import numpy as np import pyvista as pv # generate data grid for computing the values X, Y, Z = np.mgrid[-5:5:40j, -5:5:40j, -5:5:40j] values = X**2 * 0.5 + Y**2 + Z**2 * 2 # create a structured grid # (for this simple example we could've used an unstructured grid too) # note the fortran-order call to ravel()! mesh = pv.StructuredGrid(X, Y, Z) mesh.point_arrays['values'] = values.ravel(order='F') # also the active scalars # compute 3 isosurfaces isos = mesh.contour(isosurfaces=3, rng=[10, 40]) # or: mesh.contour(isosurfaces=np.linspace(10, 40, 3)) etc. # plot them interactively if you want to isos.plot(opacity=0.7) # save to stl isos.save('isosurfaces.stl')

交互绘图如下所示：

颜色对应于从标量数组中选取并由标量条指示的等值线。如果我们从文件加载回网格，我们将得到结构，但不会得到标量：

loaded = pv.read('isosurfaces.stl') loaded.plot(opacity=0.7)

缺少标量的原因是无法将数据数组导出为.stl文件：

>>> isos # original isosurface mesh PolyData (0x7fa7245a2220) N Cells: 26664 N Points: 13656 X Bounds: -4.470e+00, 4.470e+00 Y Bounds: -5.000e+00, 5.000e+00 Z Bounds: -5.000e+00, 5.000e+00 N Arrays: 3 >>> isos.point_arrays pyvista DataSetAttributes Association: POINT Contains keys: values Normals >>> isos.cell_arrays pyvista DataSetAttributes Association: CELL Contains keys: Normals >>> loaded # read back from .stl file PolyData (0x7fa7118e7d00) N Cells: 26664 N Points: 13656 X Bounds: -4.470e+00, 4.470e+00 Y Bounds: -5.000e+00, 5.000e+00 Z Bounds: -5.000e+00, 5.000e+00 N Arrays: 0

虽然原始等值面的每个都有绑定到它们的等值(提供了第一张图中看到的颜色Map)，以及点和单元法线(由于某种原因调用.save()来计算)，但在后一种情况下没有数据。尽管如此，因为您正在寻找顶点和面，所以这样做应该很好。如果需要，还可以在PyVista一侧访问这些内容，因为等值面网格是PolyData对象：

>>> isos.n_points, isos.n_cells (13656, 26664) >>> isos.points.shape # each row is a point (13656, 3) >>> isos.faces array([ 3, 0, 45, ..., 13529, 13531, 13530]) >>> isos.faces.shape (106656,)

现在，这些面孔的后勤工作有点棘手。它们都以一维整数数组的形式进行编码。在一维数组中，您总是有一个整数n告诉您给定面的大小，然后n与Points数组中的点相对应的从零开始的索引。上述等值面完全由三角形组成：

>>> isos.faces[::4] # [3 i1 i2 i3] quadruples encode faces array([3, 3, 3, ..., 3, 3, 3]) >>> isos.is_all_triangles() True

这就是为什么你会看到

>>> isos.faces.size == 4 * isos.n_cells True

你可以做isos.faces.reshape(-1, 4)得到一个二维数组，其中每一行对应一个三角面(第一列是常量3)。

【本文地址】

公司简介

联系我们